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Foreword 


Ocean  circulation  models  have  the  potential  to  provide  accurate 
forecasts,  as  well  as  nowcasts.  The  development  of  these  models  is 
facilitated  by  the  existence  of  continuous  time  series  of  quantities  to 
which  they  can  be  compared.  The  Levitus  climatology  provides  one 
standard,  but  most  useful  quantities  are  available  only  as  seasonal 
means.  Several  algorithms  are  presented  to  convert  these  into  the 
continuous  form  that  is  required. 
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Commanding  Officer 
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Executive  Summary  A  large-scale,  integrated  program  in  ocean  circulation  model  development 

has  been  pursued  by  the  Ocean  Modeling  and  Prediction  Branch,  Ocean 
Sensing  and  Prediction  Division,  Naval  Ocean  Research  and  Development 
Activity.  This  program  will  provide  the  Navy  with  accurate  nowcasts  and 
forecasts  of  the  state  of  the  ocean  on  global,  regional  and  tactical  scales. 
One  of  the  constraints  that  has  been  applied  to  this  research  is  the 
development  of  circulation  models  that  accurately  reproduce  the  measured 
climatologies  of  such  quantities  as  the  dynamic  height,  thermocline  depth, 
and  depth-averaged  density.  These  fields  are  derived  using  in  situ  data  that 
have  been  acquired  over  many  decades,  but  are  so  irregularly  distributed 
that  in  most  regions  only  annual  or  seasonal  mean  fields  have  been  derived. 
Comparisons  between  model  and  data  climatologies  are  facilitated  if  both 
are  available  as  continuous  time  series,  but  the  standard  Levitus  climatology 
is  generally  available  at  time  scales  no  shorter  than  seasonal.  This  report 
describes  several  methods  for  converting  the  existing  seasonal  quantities  into 
a  continuous  time  series  of  fields,  and  also  notes  an  amplitude  error  and 
bias  present  in  the  seasonal  values  as  presently  derived. 
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Fitting  Seasonal  Averages  with  a 
Continuous  Function 


Overview  The  Levitus  climatology  contains  annual  and  seasonal  averages  for 
potential  temperature,  salinity,  and  density  (Levitus,  1982).  To  include  this 
information  in  an  ocean  model  (either  as  a  forcing  function  or  a 
climatological  relaxation  term),  it  is  desirable  to  compute  a  more  continuous 
version  of  these  quantities.  In  the  past,  these  data  have  simply  been  inter¬ 
polated,  but  this  process  ignores  the  filtering  effects  created  by  the  seasonal 
averaging  process.  The  basic  philosophy  of  all  the  methods  described  here 
is  to  develop  a  continuous  function  whose  seasonal  averages  match  those 
given  by  the  climatology. 

The  first  part  of  this  report  describes  a  basic  problem  with  the  raw  data 
(a  bias  resulting  from  the  seasonal  averaging  and  sampling  process)  and  some 
solutions.  The  remainder  of  the  report  describes  various  methods  for 
generating  a  continuous  function  to  represent  the  climatology.  Fourier 
methods,  polynomial  interpolation,  and  splines  of  various  orders  are  analyzed 
and  compared.  The  methods  produced  similar  continuous  function,  as  well 
as  similar  first  derivatives.  The  second  derivative  in  the  case  of  the  polynomial 
and  spline  forms  is  discontinuous  between  years  or  seasons,  and  may  be 
of  concern. 

Aliasing  Present  We  can  simulate  the  process  of  developing  a  set  of  seasonal  averages  and 
in  the  Climatologies  can  show  that  there  are  two  basic  problems.  Seasonal  averaging  involves 

taking  all  the  data  that  occurs  in  a  particular  season  and  averaging  it.  The 
Levitus  climatologies  are  computed  in  this  way.  Mathematically,  this  process 
is  equivalent  to  convolving  the  original  “continuous”  time  series  with  a 
91-day  “running  average,”  then  sampling  the  resulting  smooth  time  series 
at  a  91-day  increment.  The  averaging  filter  significantly  reduces  the 
amplitudes  of  Fourier  components  whose  periods  are  shorter  than  91  days, 
but  by  no  means  entirely  removes  them.  Figures  1  and  2  show  the  results 
of  computing  seasonal  time  series  from  sine  and  cosine  waves  of  various 
periods.  When  this  averaged  time  series  is  resampled  at  a  seasonal  (91 -day) 
increment,  all  the  energy  that  remains  at  scales  shorter  than  182  days  is  aliased 
into  the  longer  scale  components.  Figures  3  and  4  show  the  response  of  the 
seasonal  averaging  filter.  All  energy  with  time  scales  shorter  than  182  days 
are  reduced  in  amplitude  by  the  91 -day  running  average.  When  these  data  are 
sampled  at  a  91 -day  increment  ,  any  energy  with  scales  shorter  than  182  days 
is  aliased  into  longer  period  components.  Once  this  aliasing  has  occurred, 
the  process  cannot  be  reversed.  The  seasonally  averaged  quantities  are 
corrupted  to  the  degree  that  significant  energy  at  scales  shorter  than  91  days 
exists  in  the  real  ocean.  The  red  spectra  of  many  oceanographic  quantities 
generally  imply  that  the  shorter  the  time  scale,  the  weaker  the  phenomenon 
(at  least  until  the  scale  is  1  day,  at  which  point  the  tidal  energy  becomes 
significant).  Exceptions  occur  at  inertial  and  tidal  periods. 

Potential  The  red  spectrum  can  be  interpreted  conversely,  as  well:  the  longer  the 
Undersampling  time  scale,  the  stronger  the  phenomenon.  Both  real,  in  situ  data  and  model 
Problems  simulations  show  that  the  sea  surface  height  (SSH)  and  the  subthermocline 
pressure  field  have  their  strongest  components  at  scales  much  longer  than 
1  year.  Figure  5  displays  the  spectrum  of  sea  surface  height  as  measured  by 
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Fitting  Seasonal  Averages  with  a  Continuous  Function 


Figure  3.  The  response  function  of  a 
91-day  running  average  (as  a  function  of 
frequency  in  units  of  1/day). 
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Figure  4.  Same  as  Figure  3,  except 
displayed  as  a  function  of  period  in  days. 
Energy  at  scales  shorter  than  91  days  is 
significant.  Any  periodicities  in  the  actual 
data  that  are  also  shorter  than  91  days  will 
be  aliased  into  longer  periods  when 
reduced  to  seasonal  means. 
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Figure  5.  Spectra  of  total  SSH  anomaly  over  a  380-day  period  comparing  IES/PC  in  situ  data  and  model-simulated  data  at  two 
locations.  Both  data  sets  show  significantly  red  character. 


an  Inverted  Echo  Sounder  (IES)  array  (Fox  et  al.,  1988)  in  the  North  Atlantic, 
compared  to  an  ocean  model  result.  Extremely  long  model  simulations  show 
that  significant  energy  exists  at  scales  as  long  at  10  years  (Fig.  6).  Any 
climatology  developed  using  only  a  few  years  of  data  will  gradually  become 
obsolete,  as  longer-term  components  treated  as  constant  begin  to  change. 
Climatologies  based  on  1-year  averages  will  be  outmoded  rapidly  since  there 
is  significant  variability  on  time  scales  longer  than  1  year.  Figure  7  displays 
climatologies  based  on  1-year  averages  of  SSH  taken  from  a  model  simula¬ 
tion  of  the  North  Atlantic  at  various  locations  compared  to  a  climatology 
based  on  all  10  years  of  the  same  data.  In  this  simulation,  interannual  varia¬ 
tions  in  the  annual  mean  sea  surface  height  are  clearly  substantial.  The  Levitus 
climatology  is  based  on  between  10  and  30  years  of  data,  but  no  data  from 
the  past  10  years  have  been  included.  If  the  variability  in  the  measured  quan¬ 
tities  is  significant  on  scales  longer  than  10  years,  then  this  climatology  may 
already  be  out  of  date  for  some  purposes. 
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Figure  6.  Spectrum  of  SSH  for  a  10-year 
model  simulation  of  the  North  Atlantic 
showing  that  the  redness  of  the  spectrum 
extends  to  scales  of  many  years. 


Biasing  of  the 
Semiannual 
Component 


After  smoothing  and  resampling  the  continuous  time  series  to  form 
seasonal  averages,  the  remaining  Fourier  components  are  the  annual  and 
semiannual  periods.  Figure  2  shows  the  effects  on  the  semiannual  component 
of  averaging  and  resampling  to  seasonal  averages.  Using  January  1  as  the 
start  of  the  time  axis  (l  =  0),  the  semiannual  sine  component  is  clearly 
reduced,  but  the  cosine  component  is  completely  filtered  out.  Even  if  the 
original  time  series  consisted  of  only  annual  and  semiannual  components, 
converting  this  data  into  seasonal  averages  destroys  the  semiannual  cosine 
wave,  so  that  the  original  simple  time  series  cannot  reconstructed.  The  sine 
component  can  be  computed,  but  not  the  cosine.  In  the  Fourier  methods 
given  here,  the  annual  and  semiannual  components  are  computed,  but  it 
should  be  emphasized  that  the  semiannual  component  is  biased  in  that  its 
phase  is  always  90°  (always  a  pure  sine  component,  with  no  cosine 
component).  The  potential  effects  are  subtle  and  not-so-subtle.  For  example, 
whatever  energy  exists  in  the  real  ocean  at  semiannual  scales  will  be  partially 
reduced,  since  only  its  sine  component  can  be  estimated.  As  an  extreme 
example,  consider  deriving  a  continuous  wind  series  to  drive  the  model  for 
seasonal  winds.  The  driving  at  semiannual  time  scales  will  be  weaker  than 
it  should  be  (on  average,  about  71%  of  what  the  total  forcing  should  be). 
An  additional  effect  of  this  bias  will  be  that  the  semiannual  forcing  will 
always  be  at  a  particular  global  phase.  A  true,  continuous  climatology  would 
almost  surely  have  a  semiannual  phase  that  varies  with  location. 

Options  for  dealing  with  this  problem  range  from  doing  nothing  to 
attempting  to  derive  the  semiannual  phase  relationship  in  some  way.  Rather 
than  include  only  the  sine  component  of  the  semiannual  climatology,  the 
semiannual  variability  could  be  totally  ignored  and  only  the  annua! 
component  used.  In  some  areas  of  the  world  (i.e.,  the  Indian  Ocean)  the 
semiannual  component  is  very  strong  and  using  only  an  annual  component 
of  the  climatology  would  significantly  reduce  the  “information  energy” 


Fitting  Seasonal  A  verages  with  a  Continuous  Function 


5 


Figure  7.  A  comparison  of  the  seasonal  climatologies  formed  from  the  IO-year-tong  data  set  of  Figure  6  and  similar  climatologies 
formed  from  individual  years  of  the  same  data  set.  Interannual  variation r  in  the  climatology  are  as  large  as  10  cm  in  this  data  set. 


being  supplied  to  the  model.  The  problem  of  transitioning  the  model 
from  being  dri  en  by  climatological  data  (to  spin  it  up)  to  driving  it  with 
real  and  current  information  will  always  be  present.  It  ir  probably  more 
important  to  include  the  most  information  possible,  regardless  of  bias,  but 
this  remains  to  be  tested  and  proven. 

Three  distinctly  different  methods  were  used  to  develop  continuous 
functions  that  represent  the  given  data  (which  consists  simply  of  four  seasonal 
averages  at  any  given  geographical  location).  In  each  case,  a  set  of  best-fit 
coefficients,  which  represent  Fourier  coefficient  weights  or  terms  in  some 
polynomial  °quation,  are  derived.  Each  grid  point  in  the  Levitus  climatology 
is  assumed  to  be  separately  treated,  which  might  resuh  in  some  undesired 
spatial  variability  (unverified).  If  this  is  the  case,  the  coefficients  could  be 
spatially  smoothed  before  being  used  to  create  the  continuous  interpolating 
functions. 

A  series  of  random  seasonal  climatologies  were  generated  r  id  fitted 
with  each  of  the  methods  described.  Three  samples  are  shown  in  Figures  8 


Overview  of 

Interpolation 

Methods 
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alues  after  being  converted  into 
• ontinuous  time  series  by  the  methods 
described  in  this  report. 
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figure  9.  The  instantaneous  time 
lerivative  of  the  curves  shown  in  Figure  8. 
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Figure  10.  Another  realization 
<see  Fig.  8). 
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Figure  11.  Another  realization 
( see  Fig.  9). 


through  13.  In  each  plot,  2  consecutive  years  of  identical  data  are  shown 
so  that  the  continuity  across  the  annual  boundary  is  evident.  The  “square 
wave”  solid  line  represents  the  seasonal  average  drawn  to  cover  its  entire 
season.  Each  example  consists  of  a  plot  of  the  data  converted  into  continuous 
form  by  each  of  the  methods  described  in  the  following  text,  followed  by 
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Figure  12.  Another  realization 
(see  Fig.  8). 
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Figure  13.  Another  realization 
( see  Fig.  9). 
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a  plot  of  the  time  derivative  of  the  continuous  interpolating  time  series.  The 
three  methods  yield  very  similar  time  series  and  similar  derivatives.  This 
similarity  is  emphasized  in  Table  1,  which  summarizes  the  results  of  35 
experiments  (similar  to  results  displayed  in  Figs.  8  through  13).  In  each  case, 
the  variability  of  the  interpolated  continuous  time  series  is  computed  for 
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Annual 

Only 

Annual  and 
Semiannual 

5th  Degree 
Polynomial 

Quadratic 

Splines 

Mean 

Variability 

STD 

Variability 

Percent 

Deviation 

0.1061 

0.1866 

0.1774 

0.1815 

0.1818 

0.0038 

2.0748 

0.3411 

0.3473 

0.3372 

0.3378 

0.3408 

0.0046 

1 .3550 

0.1851 

0.2282 

0.2292 

0.2220 

0.2265 

0.0032 

1.4091 

0.1463 

0.2057 

0.1901 

0.2001 

0.1987 

0.0065 

3.2471 

0.0697 

0.2903 

0.3203 

0.2824 

0.2977 

0.0163 

5.4776 

0.1928 

0.2128 

0.2322 

0.2069 

0.2173 

0.0108 

4.9684 

0.2352 

0.2544 

0.2768 

0.2474 

0.2595 

0.0125 

4.8280 

0.2358 

0.2472 

0.2639 

0.2404 

0.2505 

0.0099 

3.9443 

0.1856 

0.2262 

0.2411 

0.2200 

0.2291 

0.0088 

3.8490 

0.1573 

0.1984 

0.1850 

0.1930 

0.1921 

0.0055 

2.8672 

0.2385 

0.2562 

0.2665 

0.2492 

0.2573 

0.0071 

2.7656 

0.0868 

0.2166 

0.2630 

0.2107 

0.2301 

0.0234 

10.1643 

0.0257 

0.3446 

0.3857 

0.3352 

0.3552 

0.0219 

6.1734 

0.1565 

0.2350 

0.2156 

0.2286 

0.2264 

0.0081 

3.5588 

0.0937 

0.1958 

0.1963 

0.1904 

0.1942 

0.0026 

1 ,3573 

0.2023 

0.2201 

0.2324 

0.2140 

0.2222 

0.0077 

3.4493 

0.1632 

0.1931 

0.1866 

0.1878 

0.1892 

0.0028 

1.4968 

0.1862 

0.2536 

0.2611 

0.2467 

0.2538 

0.0059 

2.3254 

0.1875 

0.2769 

0.2906 

0.2694 

0.2790 

0.0088 

3.1545 

0.1180 

0.1827 

0.2109 

0.1777 

0.1905 

0.0146 

7.6534 

0.2983 

0.3149 

0.2991 

0.3063 

0.3068 

0.0065 

2.1050 

0.0649 

0.2277 

0.2386 

0.2214 

0.2292 

0.0071 

3.0901 

0.1088 

0.2394 

0.2638 

0.2328 

0.2453 

0.0133 

5.4248 

0.2701 

0.2704 

0.2747 

0.2630 

0.2694 

0.0048 

1.7839 

0.0987 

0.2404 

0.2665 

0.2338 

0.2469 

0.0141 

5.7158 

0.3058 

0.3269 

0.3316 

0.3180 

0.3255 

0.0057 

1.7361 

0.0760 

0.2453 

0.2760 

0.2386 

0.2533 

0.0163 

6.4349 

0.1113 

0.2253 

0.2726 

0.2192 

0.2390 

0.0239 

9.9833 

0.2134 

0.2236 

0.2155 

0.2175 

0.2188 

0.0035 

1.5795 

0.0559 

0.2287 

0.2384 

0.2224 

0.2299 

0.0066 

2.8630 

0.2210 

0.2278 

0.2393 

0.2215 

0.2295 

0.0074 

3.2075 

0.2065 

0.2196 

0.2158 

0.2136 

0.2163 

0.0025 

1.1481 

0.0302 

0.3235 

0.3899 

0.3147 

0.3427 

0.0336 

9.7935 

0.2471 

0.2516 

0.2456 

0.2447 

0.2473 

0.0030 

1.2290 

0.3604 

0.3701 

0.3568 

0.3600 

0.3623 

0.0057 

1.5628 

each  method  and  then  intercompared.  Excluding  the  fit  based  on  the  annual 
cycle  alone,  the  time  series  created  have  variabilities  that  deviate  from  one 
another  by  no  more  than  10%,  and  the  bulk  of  the  experiments  showed 
differences  of  only  a  few  percent. 

The  first  method  described  is  based  on  developing  a  Fourier  series  that 
represents  either  the  annual  component  alone,  or  both  the  annual  and 
semiannual  waves.  As  described,  when  the  semiannual  component  is  included 
in  the  fit,  there  is  a  bias  in  that  only  the  sine  component  can  be  estimated 
from  the  data.  However,  if  only  the  annual  component  is  used,  significant 
energy  present  in  the  semiannual  waves  is  ignored.  In  both  cases,  the  inter¬ 
polating  function  is  continuous  throughout  the  year,  as  well  as  from  year 
to  year,  and  so  are  all  its  derivatives.  This  continuity  is  important  if  the 
model  being  driven  directly  or  indirectly  by  the  climatology  responds  to 
temporal  derivatives  of  the  quantities  being  supplied.  In  the  samples  (Figs.  8 
through  13),  the  annual-only  Fit  is  shown  with  a  long  dash  line  and  the  annual 
plus  semiannual  fit  is  shown  with  a  dotted  line. 

The  second  method  is  based  on  polynomial  interpolation.  We  attempt 
to  derive  a  finite  polynomial  which,  when  averaged,  yields  the 


Table  1.  Thirty-five  realizations  of 
experiments  (see  Figs.  8-13) 
comparing  variabilities  of  the 
continuous  time  series  constructed 
by  the  various  methods  described  in 
the  text.  The  annual-only  method  is 
excluded  in  the  final  three  columns. 
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Fourier  Interpolation 


correct  seasonal  values.  The  polynomial  and  its  derivatives  are  continuous 
throughout  the  year,  but  it  is  shown  that  at  the  year-to-year  boundaries, 
continuity  constraints  must  be  relaxed  for  a  solution  to  exist.  For  example, 
in  the  case  of  a  5th  degree  polynomial,  we  can  require  continuity  of  / and 
/' ,  but  in  the  case  of  a  6th  degree  polynomial,  we  cannot  add  continuity 
of  f"  and  still  find  a  solution.  The  5th  degree  polynomial  fit  is  displayed 
in  Figures  8  through  13  using  a  dot  dashed  line. 

The  final  method  of  solution  is  based  on  splines  of  various  degrees.  Each 
season  is  represented  by  a  separate  polynomial  (linear,  quadratic,  or  cubic) 
and  then  boundary  conditions  of  continuity  are  used  to  fix  the  terms.  It 
is  shown  that  linear  and  cubic  splines  are  not  uniquely  determined  by  the 
equations  and  that  some  additional  ad  hoc  constraint  must  be  imposed. 
Quadratic  splines  can  be  derived  without  problem  (indicated  in  Figs.  8 
through  13  by  a  short  dash  line. 

In  each  of  the  methods  described  below,  the  four  known  seasonal  average 
quantities  are  represented  by  the  symbols  bx  through  bA. 

Assuming  that  the  underlying  time  series  had  a  Fourier  representation, 
the  seasonal  averaging  process  will  produce  a  time  series  that  has  a 
dc  component,  a  cosine  and  sine  component  with  an  annual  cycle,  and  a  sine 
component  with  a  semiannual  cycle.  Any  other  frequencies  are  aliased  or 
averaged  out.  The  most  general  Fourier  expansion  we  can  write  then  for 
the  original  time  series  is 

b(t)  =  a0  +  fl,cos(a>0  +  a2sin(a>f)  +  o3sin(2a>r) ,  (1) 

where  cu  is  the  frequency  corresponding  to  the  annual  cycle.  If  time  is 
measured  in  days,  then  the  first  day  (Jan  1)  corresponds  to  t  -  0  and  the 
fundamental  frequency  will  be  given  by  co  =  2tt/364  (for  the  case  of  a  model 
year  with  364  days).  If  time  is  measured  in  weeks,  then  co  =  2tt/52  and 
the  first  week  will  be  at  t  =  0. 

The  four  seasonal  averages  of  Equation  (1)  are  given  by: 


bl  =  <  b(t )  >,  =  a0  +  a]  <  cos(a >/)  >, 

-f -  .  .  .  +  a3  <  sin(2c ot)  >,  ,  (2) 

b2  =  <  b(t)  >2  =  a0  +  <  cos(a U)  >2 

+  .  .  .  +  a3  <  sin(2 wt)  >  2 ,  (3) 

b3  =  <  b{t)  >3  =  a0  +  a,  <  cos  (tot)  >3 

+  .  .  .  +  o3  <  sin(2a)f)  >  3 ,  (4) 

bA  =  <  b(t)  >4  =  a0  +  c,  <  cos(cof)  >4 

+  .  .  .  +  a3  <  sin(2cof)  >4  ,  (5) 


where  6,  represents  the  January-March  average,  b2  the  April-June  average, 
etc.,  and  the  symbol  <...>,  means  the  average  for  season  i. 
Performing  the  various  averages  as  continuous  integrations  of  the  sinusoids 
provides  the  key  to  the  simplicity  of  the  final  algorithm:  all  the  various 
averages  reduce  to  either  +  2/n  or  -2/n. 
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The  fina'  seasonal  average  equations  are,  simply, 


hj  =  a0  +  (2/n)  (+  a,  +  a2  +  a3)  , 

(6) 

b2  =  a0  +  (2/n)  (-  a,  +  a2-  fl3)  , 

(7) 

b2  =  a0  +  (2/n)  (-  a,  -  a2  +  a3)  , 

(8) 

bA  =  a0  +  (2/n)  (  +  a,  -  a2  -  a3)  . 

(9) 

Adding  all  four  equation0  Produces  the  expected  value  for  the  constant 
term. 

ao  =  (^i  +  b2  +  bi  +  bA)  . 

(10) 

Adding  the  resulting  equations  in  various  pairs  produces  the  final  estimates 
for  the  remaining  a’s: 

ai  =  (rt/8)  (Z>,  -  b2-  b3  +  bA) , 

(11) 

a2  =  (n/8)  (6,  +  b2  -  b3  -  bA) , 

(12) 

a}  =  (n/8)  ( b ,  -  b2  +  b3  -  b4 )  . 

(13) 

Thus,  from  the  original  seasonal  amplitudes  bjt  we  can  immediately 
construct  the  coefficients  a,  in  the  expansion.  Seasonally  averaging  the  time 
series  will  reproduce  the  given  b ,  exactly. 

If  we  choose  to  ignore  the  semiannual  component,  we  can  use  the  fact 
that  this  component  is  orthogonal  to  the  annual  one,  so  that  simply  setting 
a3  =  0  will  give  the  least-squares  best  fit.  Alternately,  we  can  go  through 
the  mechanics  and  verify  that  the  constant  and  annual  cosine  and  sine  terms 
have  the  same  coefficients  as  derived  above. 

One  obvious  way  to  interpolate  the  given  seasonal  averages  is  to  use  a  Polynomial  Interpolation 
finite  polynomial.  That  is,  we  attempt  to  find  a  polynomial  of  degree  N 
such  that  the  time  series 

b(t)  =  I  c,t‘  (14) 

<  =  o 

integrated  over  3 -month  periods  provides  the  proper  seasonal  averages.  The 
additional  constraint  in  the  problem,  however,  is  that  the  time  series  must 
have  a  period  of  1  year,  and  it  should  be  continuous  (and  have  continuous 
derivatives)  everywhere,  including  where  the  function  shifts  from  the  end 
of  the  year  to  the  beginning  of  the  next  year.  If  we  normalize  the  time  variable 


to  the  range  [0,1], 

then  we  have  the  additional  set  of  N  constraints: 

b(  0)  = 

=  b(  1) , 

(15) 

b'(  0) 

=  fr'(l). 

(16) 

b"(0) 

=  b"(  1), 

(17) 

•  •  •  » 

(18) 
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b "(0)  =  b'N-  »(1) . 


(19) 


Polynomial  Interpolation 
with  Relaxed  Constraints 


A  final  constraint, 

b<»>(0)  =  b<">(l)  ,  (20) 

is  satisfied  automatically  because,  it  simply  requires  that  cN  =  cN,  but 
problems  are  created.  We  have  N  +  1  unknowns  (c,),  but  N  +  4 
constraints  (the  N  equations  in  (20),  plus  the  four  seasonal  averages  that 
are  required).  One  might  be  tempted  to  solve  this  as  a  least-squares  problem, 
but  first  it  will  be  shown  that  the  derivative  constraints  prevent  a  solution 
from  being  found. 

As  an  example,  take  the  simple  quadratic  function, 

/(/)  =  a  +  bt  +  ct2  .  (21) 

The  time  variable  is  normalized  to  go  from  t  =  0  to  t  -  1,  so  the 
requirement  of  continuity  of  the  function  and  its  derivatives  implies 
the  following  set  of  equations  (note  first  that  /'(/)  =  b  +  2c/): 

m  =  /(l)  -  a  =  a  +  b  +  c ,  (22) 

/'(0)  =/'(l)-*  b  =  b  +  2c.  (23) 

The  last  equation  sets  c  to  zero.  The  previous  equation  then  requires  that 
b  be  zero  as  well.  Generally,  one  can  show  that  the  only  finite  polynomial 
function,  which  is  periodic  and  continuous  with  continuous  derivatives  that 
match  across  the  end  points,  is  the  function ./(/)  =  constant.  (If  we  want 
an  infinite  polynomial,  all  we  have  to  do  is  expand  the  Fourier  series  derived 
above  for  the  annual  and  semiannual  fit.) 

As  seen  in  the  last  section,  no  single,  finite  polynomial  will  interpolate 
the  seasonal  data  in  a  satisfactory  way.  The  constraint  of  having  all  its 
derivatives  continuous  across  yearly  boundaries  would  have  to  be  removed 
to  permit  a  solution.  A  discontinuous  derivative  in  the  time  domain  implies 
a  significant  amount  of  power  at  high  frequencies  in  the  Fourier  domain, 
which  is  undesirable,  but  we  can  still  generate  the  polynomial  and  assess 
its  accuracy. 

We  have  four  constraint  equations  that  come  from  the  given  seasonal 
averages,  and  we  probably  want  the  function  and  at  least  its  first  derivative 
to  be  continuous  across  yearly  boundaries.  This  operation  provides  six 
constraints,  so  the  simplest  polynomial  will  have  the  form, 

5 

b(t)  =  I  cV .  (24) 

1  =  0 

It  is  convenient  to  define  the  time  to  be  in  the  range  [-2,2].  That  is,  /  =  -2 
represents  the  beginning  of  January  1  and  t  =  +  2  would  be  the  end  of 
December  31 .  Integrating  b(t)  from  /  =  -2  to  /  =  -1  yields  the  average  for 
the  first  3  months  of  the  year,  which  should  equal  b integrating  from 
/  =  -1  to  /  =  0  should  yield  the  given  average  for  the  second  3  months 
of  the  year  b2,  and  so  on.  Continuity  of  b(t)  and  b'(t)  yields  the  final  two 
constraints.  These  six  equations  in  the  six  unknowns  are  represented  by  the 
matrix  equation. 
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fl 

-3/2 

7/3 

-15/4 

31/5 

-21/2 

127/7  > 

ro 

fM 

1 

-1/2 

1/3 

-1/4 

1/5 

-1/6 

1/7 

c. 

b2 

1 

1/2 

1/3 

1/4 

1/5 

1/6 

1/7 

C2 

b, 

1 

3/2 

7/3 

15/4 

31/5 

21/2 

127/7 

C3 

— 

b4 

0 

1 

0 

4 

0 

16 

0 

0 

1° 

0 

1 

0 

8 

0 

48  J 

lC5  j 

1°, 

(25) 


(25) 

This  matrix  can  be  inverted  exactly  to  yield: 

Cq  =  (-37/300)0,  +  b4)  +  (187/300)(b,  +  b}) ,  (26) 

c,  =  (5/18)0,  -  b4)  -  (3/2 )(b2  -  b2)  ,  (27) 

c2  =  (2/5)0,  -  b2-  b2  +  b4)  ,  (28) 

c3  =  (-47/72)0,  -  b4 )  +  (9/8)02  -  b2)  ,  (29) 

c4  =  (-1/20)0,  -  b2  -  b,  +  b4)  ,  (30) 

c5  =  (7/48)0,  -  b4)  -  (3/16)02  -  b})  .  (31) 


At  each  grid  point  in  the  model,  then,  the  four  seasonal  average  quantities 
are  transformed  into  a  set  of  six  coefficients  for  a  5th  degree  interpolating 
polynomial. 

In  the  normal  method  of  spline  interpolation,  we  are  given  a  set  of  data 
points,  and  we  desire  a  set  of  cubic  polynomials  defined  between  the  points 
that  match  at  the  points,  as  well  as  having  matching  first  and  second 
derivatives  at  the  points.  In  our  problem,  the  points  represent  averages  over 
a  season;  thus,  we  are  more  interested  in  finding  a  set  of  interpolating 
polynomials  that  match  well  at  the  boundaries,  and  that  yield  the  proper 
seasonal  averages  when  appropriately  integrated.  This  philosophy  is  similar 
to  the  usual  definition,  so  the  method  will  still  be  referred  to  as  spline 
interpolation. 

The  first  spline  method  discussed  will  simply  use  linear  functions  in  each 
season.  That  is,  in  each  quarter,  the  function  is  assumed  to  have  the  form 
bj(t)  =  c0,  +  cut.  This  gives  eight  unknowns.  The  seasonal  averages 
provide  four  constraint  equations,  and  the  requirement  that  the  various 
functions  match  at  their  respective  boundaries  provides  another  four.  These 
constraints  would  appear  to  provide  equations  to  determine  the  unknown, 
but  the  equations  are  not  independent. 

Let  the  time  argument  t  vary  from  -1  to  +  1  in  each  season  (quarter). 
The  seasonal  averaging  process  (integrating  each  function  over  [-1,1]  yields 
the  constant  terms  directly. 

c0(  =  0.5  b , .  (32) 

The  requirement  of  continuity  across  seasonal  boundaries  (including 
between  seasons  4  and  1  at  the  year-to-year  boundary)  produces  (after 
rearrangement  and  replacing  c0(). 

c„  +  cu  =  0.502  -  &,) ,  (33) 


Spline  Interpolation: 
Introduction 


Linear  Splines 
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(34) 


Cubic  Splines 


Quadratic  Splines 


Cj2  +  Cjj  —  0.5(63  6^)  • 

C13  +  C14  =  °-5(64  -  .  (35) 

Cu  +  c„  =  0.5(6,  -  64)  .  (36) 

Unfortunately,  these  equations  are  inconsistent,  as  can  be  demonstrated 
by  adding  the  first  and  third  equations  and  subtracting  the  second  and  fourth 
equations.  This  leads  to  an  unacceptable  constraint  on  the  data  itself;  that  is: 

0  =  6,  +  63  -  62  -  64  .  (37) 

We  can  solve  problems  of  this  type  if  we  bring  in  an  additional  ad  hoc 
requirement.  For  example,  we  may  try  to  find  the  set  of  linear  functions 
that  has  the  least  total  variance  while  still  exactly  satisfying  a  subset  of  the 
original  constraint  equations,  or  we  might  look  for  the  set  that  has 
the  smallest  overall  coefficients  in  the  fit.  This,  then,  becomes  a  constrained 
least-squares  problem  and  methods  of  solving  them  (particularly  with  the 
small  numbers  of  equations  involved)  are  no  more  difficult  than  normal 
least  squares  problems  (Brandt,  1970;  Claerbout,  1976;  Golub,  1983; 
Meyer,  1982).  In  particular,  if  the  absolute  constraints  we  require  are  given 
by  the  matrix  equation 

Gc  =  d  (38) 

where  c  is  the  column  of  coefficients  we  are  looking  for,  and  the  goal  is 
to  minimize  the  energy  of  the  coefficients  (that  is  minimize  the  scalar  cTc), 
then  the  solution  is  given  by: 

c  =  GT(GGT)-'d  .  (39) 


As  above,  we  attempt  to  find  a  set  of  four  cubic  polynomials,  one  for 
each  season,  which  match  at  the  seasonal  and  annual  boundaries  and  also 
have  continuous  first  and  second  derivatives  there.  Further,  we  require  that 
the  seasonal  averages  computed  from  these  functions  match  the  four  given 
values.  As  above,  we  start  off  appearing  to  have  enough  equations  to  solve  for 
all  16  unknowns  (four  seasons,  and  four  coefficients  for  each  cubic 
polynomial):  the  four  seasonal  average  equations  and  four  equations  each 
from  continuity  of  the  functions,  their  first  and  their  second  derivatives  across 
boundaries. 

In  fact,  again  as  above,  the  equations  turn  out  to  be  linearly  dependent, 
and  some  type  of  least-squares  solution,  constrained  to  meet  the  remaining 
independent  equations  exactly,  must  be  performed.  Rather  than  do  this 
solution  at  this  time,  the  requirement  of  continuity  of  the  second  derivative 
is  dropped,  and  the  order  of  the  interpolating  polynomials  is  dropped  down 
to  quadratic  in  the  following  section. 

After  reading  the  previous  two  sections,  one  might  wonder  why  another 
spline  order  should  be  attempted,  but  in  this  case  it  can  be  solved  without 
external  ad  hoc  requirements  or  least-squares  criteria. 

Each  of  the  four  seasons  is  represented  by  a  unique  quadratic  function 
in  t.  For  convenience,  separately  define  the  time  to  run  from  t  -  -1  to 
t  =  +  1  for  each  season.  Then  the  four  functions  each  have  the  following 
form 
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/,(')  =  %  +  Cut  +  c2,t2  ,  (40) 

flit)  =  cu  +  2c2/t  .  (41) 

The  continuity  of  /  across  seasonal  boundaries  provides  the  first  four 
equations. 

/](+  1)  =  ~ *  C01  +  C11  +  C2\  =  C02  ~  C12  +  C22  >  (42) 

•/■>(  +  0  =  /3(-l)  *  c02  C12  ^22  =  C03  ~  C13  ^23  ’  (43) 

/3(+  D  =  /«H)  —  c03  +  c13  +  c23  =  -  c14  +  c24  ,  (44) 

y"4(  +  1)  =  /i(~l)  ^04  "*■  C14  C24  =  C01  ~  CI1  C2l  ’  (45) 

Continuity  of  the  first  temporal  derivative /'  across  seasonal  boundaries 
implies  that 

/,'( +  1)  =  /i(-l)  -  c„  +  2c21  =  c12  -  2c22  ,  (46) 

/2'(  +  1)  =  /3(-l)  -  ci2  +  2cn  =  c13  -  2c23  ,  (47) 

/3'(+  1)  =  /J(-l)  -  c13  +  2c23  =  c,4  -  2ClA  ,  (48) 

/i(+  1)  =  /,'(-D  -  c14  +  2cm  =  c„  -  2c21  .  (49) 

The  final  four  relationships  among  the  12  c,2  are  derived  using  the 
seasonal  averages.  For  each  of  the  four  seasons,  there  is  a  constraint  that 

J  ft(t)dt  =  b,  ,  (50) 


where  the  bi  are  the  four  given  seasonal  means.  Performing  the  integration 
results  in  the  four  equations: 

3co,  +  c2/  =  36<  •  (51> 


All  of  these  relations  can  be  compactly  represented  in  matrix  form: 
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This  matrix  can  be  solved  to  yield 


19 ,  1 


1  .  1 


Loi 


c.,  = 


«(b2-b4) 


C21  =  - 


9  .  3  3  3  , 

T66>  +  I^-T6^  +  I^ 


(53) 

(54) 

(55) 


for  the  first  season.  By  direct  solution,  or  by  symmetry,  the  solutions  for 
the  other  three  seasons  can  be  found  by  permuting  the  index  to  represent  the 
season  appropriately.  For  example,  to  solve  for  season  2,  replace  the  second 
index  of  the  ctJ  by  “2”  and  permute  the  subscripts  on  b  from  (1234j  to 
[4123].  By  continuing  to  permute  circularly,  the  values  of  c/2  through  c/4 
are  generated. 

Summary  Several  methods  permit  quantities  known  only  as  seasonal  values  to  be 
converted  into  continuous  time  series.  These  fields  can  be  used  as  continuous 
forcing  functions  for  models  or  as  a  continuous  climatological  relaxation 
term  in  the  model  equations.  Errors  in  the  seasonal  values  due  to  averaging 
and  sampling  (aliasing)  were  described.  A  companion  report  is  being  written 
that  proposes  a  method  for  computing  continuous  time  series  that  do  not 
have  these  undesirable  effects. 
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